


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1971 


Numerical analysis of multigroup neutron flux 
in a bare fast reactor 


Bosworth, Robin. 


Monterey, California ; Naval Postgraduate School 


http://ndl.handle.net/10945/15831 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
D U DLEY research materials and institutional publications created by the NPS community. 
get Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS'‘s first 
KNOX appointed — and published — scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 


NUMERICAL ANALYSIS 


OF MULTIGROUP NEUTRON FLUX 
IN A BARE FAST REACTOR 


Robin Bosworth 

















Schoo! 





Ye F.4 west e LE Mae: 
THESIS 


| { Sere we Gyr er er) re ee is oe « - Sau - ao ee ee OSE 2 Beem 








NUMERICAL ANALYSIS 
OF 
MUERIeROUr NEUTRON TU Mine Ban Ast REACTOR 






by 








Robin Bosworth 


Thesis Advisor: 





Approved for public release; distribution untlimcted. 


Q"y 


- 
see a oe me eer me ce 








Numerical Analysis 
of 


Multigroup Neutron Flux in a Bare Fast Reactor 


by 


Robin Bosworth 
Lieutenant, United States Navy 
B.S., United States Naval Academy, 1964 


Submitted in partial fulfillmert of the 
requirements for the degree of 


MECHANICAL ENGINEER 


' from the 


NAV eros! CkADpuUst EE SBOnOOL 
June 1971 


es 
7 26.7 
ost 





ABSTRACT 


This study consists of the development of a computer program to 
numerically solve the space and energy dependent multigroup neutron 
diffusion equations in a bare homogeneous fast reactor core or reactor 
material assembly. 

The resulting program is unique in that it was designed for future 
use by Naval Postgraduate School students undertaking experimental 
studies in neutron diffusion with limited time to determine numerical 
solutions for verification of their results. 

The equations are solved iteratively in cylindrical geometry using 
a point successive overrelaxation technique. Convergence between 
successive iterations was less than 10° after fifty iterations. 

The program was tested using ANL three group data. Flux shapes 
and energy spectra were determined for a typical fast reactor core 
and for a solid iron cylinder with a source at its center. ‘the program 
was also used to determine criticality. 

Computation times were from one to ten minutes with less than 


150K words of core storage using the IBM 360/67. 
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PeINTRODUC Hey 


A. BACKGROUND 

It has been found that the energy-dependent diffusion equation 
provides a fairly accurate description of the neutron behavior in fast 
reactor systems [1]. 

One way of conveniently using this equation is the now familiar 
multigroup method which was developed for reactor calculations by 
F., L. Friedman, A. M. Weinberg, and J. A. Wheeler [2] about 1947. 
Since that time, many solutions to the age-diffusion equation have 
resulted from the use of the multigroup method. Some analytical 
and semianalytical solutions have been carried out for a few group 
equations [3] but most of the work done has been numerical due mainly 
to the need for large numbers of energy groups, spatial points, and 
spatial regions. 

The numerical solutions to unique problems and classes ot 
problems have evolved in the form of computer programs or codes. 
The degree of complexity and sophistication of these programs has 
roughly paralleled that of the digital computer. At this point in time, 
the literature abounds with numerous and varied codes and references 
thereto. Sangren [4] states that by 1959 there existed in.the United 
States more than 300 nuclear reactor codes for digital computers. 
Various groups have classified and categorized these codes [4]. In 
1961 the American Nuclear Society and the Argonne National Laboratory 
established a code center for the purpose of collecting and disseminat- 
ing the multitude of codes that existed and for the future codes which 
were sure to be forthcoming [5]. Around 1965 the European counter- 
part to this center was established [5]. 

Most of the more recent codes (those published after 1965) 


examined in connection with this study are somewhat specialized 








along the following lines: a) number of energy (lethargy) groups; 

b) number of spatial dimensions; c) number of regions; d) cross-section 
data input: e) geometry; f) number of mesh points; g) method of equation 
solution; h) finite difference approximation and i) many other features 


related to poisons, burn-up, perturbations, computer system, and others. 


Bee OBJECTIVES 

This study was undertaken to provide a relatively simple-to-apply 
multigroup diffusion program for application to fast reactor cores and 
materials. in the future students at the Naval Postgraduate School 
will be able to use this program in conjunction with experimental 
neutron diffusion studies. | In view of the large number of complex 
codes that exist ‘and the difficulties related to modifying these codes 
for use on a specific problem and a specific computer system, itis 
anticipated that valuable research time may be saved by use of this 
Prcoram. 

Once the program was developed, initial testing would be done 
using a subcritical assembly. The first phase of this test would 
involve the solution of the three group diffusion equation in a typical 
fast core composition. From this, the energy and spatial dependence 
of the flux could be determined and comparisons made with the solution 
to the diffusion problem wherein space and energy dependence are 
treated as separable variables. 

The second phase of the initial test would involve the problem 
of the typical core assembly with a Pu-Be neutron source located at 
its center. This would then be modified to a solid iron cylinder with 
a centrally located Pu-Be neutron source. The flux shapes and 
energy spectra resulting from this study could then be used for 
comparison with experimental data resulting from local research. 


Experiments in this area are currently underway with applications 








to fast reactor cross section data evaluation, materials and shielding 
ReokenenGe soi me 

Using the program to determine a critical composition and size for 
the typical fast reactor was another desired goal. In this area, a 
study was to be made on the effect mesh spacing has on criticality 
Hany . 

Other objectives were to compare the energy spectra of the three 
and the sixteen group solutions and to ascertain any difference in 
the flux shapes and spectra by using the ANL sixteen group data 
and the LASL sixteen group data. 
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A. THE BOLTZMANN TRANSPORT EQUATION 

The conservation or balance of neutrons in a nuclear reactor can 
be completely described only by simultaneously specifying the distri- 
bution in space, time, energy and direction of motion of the neutrons. 
The general equation governing this balance is called the Boltzmann 


transport equation and may be written as 


PRODUCTION - LEAKAGE - ABSORPTION = 


Ae 
d|> 


an 


where Se: 


is the time rate of change of neutron density. 


The neutron density is in general described as 


number of neutrons at time t in volume 

dr about r, whose speeds lie in dv about 

v and whose directions of notion (velccily 
vectors) lie in the differential solid angle 
| angle do about direction fA. . 


HH 
Fr 


and neutron flux is defined as 
PCF, nr, HL, t= wnci, nv, A,t) 1-3 


Meghreblian and Holmes 1 describe Equation I]-1 verbally as 


number of neutrons at speed in dv about 
I |v moving in direction °% which appear 
per unit time from source in dr 


number of neutrons of direction © gained 
per UNLE time in Gr from Scattering 
+I] | collisions which scatter the neutron from 
- lall directions 7’ to -% and all speeds 
to speed in dv about v 


dtc 








net number of neutrons moving in direction 
~ at speed in dv about v that are lost 
per unit time by leakage through the 
boundaries of dr : 


ll 


moving in direction — which are removed 
per unit time from dr by absorption 
collisions 


= 


Fe of neutrons at speed in dv about 


moving in direction ™ which are remove 
per unit time from dr by being scattered 


number of neutrons at speed in dv about 
d 
into a new direction 


=VI 


change per unit time in net number of 
neutrons of direction - indr 


The first term in Equation IJ-4 accounts for the rate at which 
Memions with speed in dv about v ane with direction in dO. annear 


from sources in dr. This term is then written as 
5 Crow, ans, 2.) cimeimace@aan II-5 


The second term of Equation II-4 gives the total number of neutrons 
which appear in the differential element drdvd per unit re 
with speed in dv about v and with directiondR about M by 
being scattered from all possible speeds v and direction of motion 
The number of neutrons in drdvd FAL which expemenqe tis scarier 


is 


Elm DR, vj) Bt) dFdw'd A’ | 1-6 








The probability that the scattering collision results in the neutron 
having speed in dv about v and directiond ® about M™ is 


defined by 


CMa Os Ne el) = 
probability that a neutron with an ) 
initial speed v and direction of 
motion ©" , when scattered, emerges 
from the collision with a speed in 
dv about v and direction of motion 
indm about fi 


ay 


If Equations IJ-6 and JI-7 are multiplied together and integrated over 
all possible v' and ’ , the result is the total gain of neutrons in 
the differential element dfFdvdf_ due to the scattering in of neutrons 


from higher speeds and different directions of motion written 


PCF, a a0! cep) reels eC fie ae 
Gin Gt GG) ceo) aaa 


so 


The third term of Equation JI-4 represents the net loss of neutrons from 
drdt due to leakage through the elemental boundaries of dr and is 


equal to 
IL-VPdFdA 11-9 


A rigorous derivation of Equation II-9 is given in reference [1]. 
Term number four of Equation I]-4 is the absorption term and accounts 
for all neutrons lost from drdvd 1 per unit time by absorption; it 


is given by 


Faw) PCF ny, H,t)dFdwd A ; I-10 


Ts 








The rate at which neutrons are scattered out of drdvdit into new 
directions fi’ is the fifth term of Equation II-4 and is written as 


Zs Cw) PCF mw, Bt ydFdwda . I-11 
The final term of Equation II-4 is the time-rate of change of the 
neutron density in drdvd This term can be written as 
on _ dee = Lee 
Spo SA Al GE.) Valera: eae = “ein 
1 oO — — _ = a 
rape OCT hy LOU te) GIT > cl ies 


By placing Equations II-5 through JI-12 into the verbal Equation II-4 
the result is the Boltzmann 


and dividing out the differential drdvd 
transport equation for the distribution of neutrons of all energies in 


Space and time 
at #¢ r,t) +2 .(w) PCF, Ht) 
ans pelea) = SC, Core) 

Car eOmnc VER) SIC, Mai, oe elie 60." 


an 


Tha ths 


z.. me 


The second term in Equation II-13 is the total loss of neutrons due to 
The total 


scattering out to new spaces and directions plus absorption 


eve@ss section ban Cv) is used. 
Equation II-13 is perfectly general with the following two 
) the medium is assumed to be isotropic and homogeneous 


exceptions: 
and b) neutron slowing down is due only to elastic-scattering collisions 


can nuclei. 
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of the gth group is defined by the integral 


Ee 4-1 
Do CF) - | OCF,E) dE II-17 
Eg E 
where the upper limit on the integral E q-| is the lower boundary 


of the energy interval and the lower limit on the integral Eg is the 
upper boundary of the energy integral. The energy-dependent flux at 
the point F is Olea e) 

Within each group the diffusion of neutrons is described by the 


average diffusion coefficient. For simple diffusion theory this is 


written as 
\ 
De = _ II-18 
S S tr 9 
wnere tr, is the group macroscopic transnart cross section and 
may be obtained from 
= \Z < I-19 
De: outa, g 2 N Ort r,3 ) 


where the summation is over all elemental components and the group 


a 8 8 Zz. t 
microscopic cross section OLY 9 in the absence of pronounced 
d 


resonance effects due to capture and/or.scattering with the energy 


group is given by 


EgH 
Zz 
ae Fes oe Fen les) DCe) dE 
a I-20 


— 
i 3" BCE) dE 
Ei 








B. THE AGE-DIFFUSION APPROXIMATION 

The age-diffusion approximation to Equation II-13 is developed 
in reference [1] by expanding the Boltzmann transport equation in 
Spmecricall harmonies. 

In the steady state the resulting first order approximation to 


Bouation I]=~13 is 


-DIW VE DCRU) + EBaq tu) DCF,W) = SFU) 


- Ji he 
a ae a ee 


Since the change in the slowing down density of neutrons within the lethargy 
interval du must be balanced by the flow of neutrons scattered into 


du from elsewhere, the last term of Equation IIJ-14 may be written as 


S ens a 
7 BE UCU) = | Zs (UU) BU) Av 


And Equation JI-14 may be rewritten substituting energy for lethargy 


variable as 
~-DCE) U° DCH ED + SCE) PF,E)D = 
I-16 


S CR) eas ‘ He OS ey 2) Ce cis” 


Complete details of this analysis may be found in reference [1]. 


C. MULTIGROUP FORMULATION 
In forming the multigroup diffusion equation, the initial step is 
to divide the complete neutron energy (lethargy) range into N groups 


which may or may not be of the same width. The neutron flux 


1€ 








Neutrons may be removed from a particular group g through an 


absorption interaction described by the group absorption cross section 


Fag 
2654) > 2 6.q * 2 cae ee ee ) It-21 
where 
Pane = fission cross section of group g 
2 Byls = capture cross section of group g 
Zi mg = > Pe nsacae = inelastic removal cross 1-29 
We section from group g 
ms erg = 2 ane elastic removal cross 
K#9Q section from group g 


The group macroscopic cross sections defined in Equation II-2] above 
are obtained irom the corresponding group microscopic cross sections 


for each distinct element or isotope in the reactor such that 


2414 “2, N Tyg > II-23 


Z 
where q is the particular reaction and N’ equals the number of atoms 
of material z per cubic centimeter. 

The multigroup microscopic cross sections for each material are 


defined in reference [7] as follows: 


GROUP CAPTURE 


Eon 
[ [ot (E)+ oY 2 (Et J D(EVdE 
Zz. Egt 24 


Fosqg = Eat 
a {° OC EGE 
Ege 
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GROUP FISSION 
Egt - 
Tris (E) PCE)dE 


E 
| " GCENGE 
Egt 


Fa E gt 
S59 = y 11-25 


GROUP SCATTERING (IN OR OUT) 
Eot Egi z 
[ a2 CE FE) BCE) dE dE’ 


. E'sEgitg-egt inn 
a. = 
IN = Eay 
et 3 K it S P(E) ae 
Eat 


where ae and on are the upper and lower limits of group g respectively 
and D (E) is the flux per unit energy interval at E. 


Neutrons may be added to a particular group g by a fission source 


N 
S69 = Xo > (M Zsdg Do : 11-27 
: Q=I 


where X g is the fraction of fission neutrons born into group g such that 


Egtt 
Xo = PECayCle 5 1-28 
Boe 
(Eee = 2. Ne (DIG 5 I1-29 
and oie ie > 
[ a? ots Ce) Pe) de. 
ao Me gt 
(Dog) = Eqn II-30 
[-* dende 


Bae 








The steady state diffusion equation for the gth energy group, where 


g represents any group between one and N, can now be written 


_ Ds VAs (FI ~ Zags Py) 
[2 Eg +hpah)+ ao. ZTCh+9) Dy (F) 
u X > hE sn Dr (F) Sg CF ) ccs” 


In Equation II-31 the inelastic removal cross section from group g 


appears as the third term of the equation rather than as an implied 

part of the second term. This equation now represents a set of N 
coupled partial differential equations independent of time and applicable 
to one region of areactor. The right hand side of the equation 
represents the number of neutrons at energies within group g coming 


from any extraneous neutron sources. 


D. SPACE-INDEPENDENT FORM OF MULTIGROUP EQUATION 

For this case, it is assumed that the extrapolation distance at the 
surface of the reactor is identical for all energy groups and thus is 
independent of energy. If this is true then all energy groups have the 


same spatial dependence. The flux is then written as 
Dg(F) = Dg FCF) ree 


where Pg is the magnitude of the flux in the gth energy group, and 


the function F(r) satisfies the equation 
Pi Pe = 
V ar) > BIaCe) =O MSOs 


the factor Be is the square of the first eigenvalue of Equation II-33 


and is called the buckling. 
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By substituting Equations II-32 and 11-33 into Equation [I-31] 
expanding the resulting equation in a series of eigenfunctions, and 
applying the orthogonality condition, the result will be a set of N 


linear algebraic equations in N unknowns By ; do, -*- By written 


“EDs B° + 2ayq SS = Cah) | Pe 
DOE tng) Bn +g Dn Ean Gr > = i 


fie 


This set of equations is easily solved by Cramer's Rule. 


pee EXPLICIT SPATIAL DEPENDENCE 

In order to solve the set of Equations II-31 by numerical methods, 
the Laplacian operator in the first term must be replaced by its equiva-: 
lent finite difference approximation. The resulting equation for a 


typical energy group g in cylindrical geometry is 

Dg kt ; De (r+i, 2,6)-2Qq (r,z,0)+Bg(r-i,z,e) $ 
tom § Bg (r+i,z,0)~ So (r-1,2,6)3 +72 $ Po Cr, Z-1,0) 

—2 9 (¥,z,e) + Palr z+ 03+, § Pa(nz, O+!) 
~ 299 Cr, Z,O)+ Bg (r,Z,0- -)3|- 5 7-349 Ba Cay Op 
ae (g+h) | DgqCr,Z,6) + > Z2.Ch>q) Bn(r,z,6) “°° 
+X 20h Zsyn Pry, Z ae =Sq(r,z, ©) 


Appendix A contains the details of the derivation of the finite difference 


approximation and the definition of the symbols used in Equation II-35. 
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Throughout this investigation, circular summetry was imposed. 
Thus the third term in Equation II-35 is zero. For mesh points not 
on the longitudinal axis (See Appendix A for finite difference operator 


used on this axis.), Equation II-35 becomes 


Dg [a2 { Bg Crt, 2-249 Cr 2) + 3 CaN 72) @ 

torn f dg (r+, 2) - Dar, 2% +a F Yer, z+ 
“digs Z)+ fg Cr, z-1)3 | - Zo, 9g fa (2) oat 
a 2 ZCg+h)| pg cr,z) + = th+9) Sncr,z) 


‘ > 2) Pr a 1 @ anti 2) SST.“ GO 20> 
n=l 


F, MATRIX FORM OF MULTIGROUP EQUATIONS 

For the purpose of numerical analysis and computational ease, 
Equation II-36 is reduced to the matrix form. In Appendix B, this 
equation is first expanded to fit the available sixteen group data, 
then reduced to a form more suitable for matrix notation. Finally, 


it is reduced to the form | 
Al Bg (r~1, Z)4 A2 Gg (r,Z) AB Pa Cr 41,2) 
Wate OG, 2-1) +S Ga Ca 241) = Saath, zZ) 
EAS, D, Ce, Z ee AS gaia 


+ AFg Gao (r,Z) +-°-- +AFis Pic (rz) 
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fo reduce this set of multigroup equations to matrix form, the 


spatial dependence is considered first. 


Ie rtixeviecEor do of group g is defined as 


Pg (151) 


Bg (L, 1) 


s 


ie 


Si 


Pg (1, M) 


Bo(Lm) | 


where @g(L,2+m) and Gg (2~+L,mM) are on the extrapolated 
boundary and are set equal to zero. The fluxes GC\,X) and ZBCX,}) 
owe set equal to the fluxes @(€3,X) and BCX,3) respectively 
in order to satisfy the normal derivative boundary condition. See 
Appendix C for mesh point numbering scheme. 


The matrix of coefficients Ag is written 


OFA O27 @ AlA2 ise] 6 ee 5076 76 
OO A4:--O OAIAZA3::O O CASO:::O 


5 II-38 


The terms on the right hand side of equation B-4 are called the 


source terms and Wa is the source vector for group g. It can be 


(be 








represented in matrix form by 


-> Sy Br 2 AF, Bh + Sq Sh 


h=Gr 


where the flux vector Diy is given as 


OGD 


Pr Ca) 


Pr i Oy (L, 1) : 


Bh (L,M) 


the scattering matrix ASh 


pan 5@ OW NCIN@r © .Omineee 
OO oy oO a 


Ii-41 
ES. 


= ee ee ee oe 


the fission matrix AFh 


OTe. OO “ae hae. @ oo 

OO: 16" One eS Ce = tae 
AF, = 

O O OF: O O 
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I1-40 


ae 
1S) 
ti=42 


AF, | 








and the extraneous source matrix Sg 


Sq = - 11-43 
ee OS ee ee cS, OP =. 
=| 


The result is then a set of equations, one for each spatial mesh 


point for the energy group g. These may be written as 


Ag Da =Wa | 11-44 


Next the equations are cast according to the group dependence. 
Define gd as 


pb i 11-45 


PY Ome ‘oe 
O Agr : O 

A ae > I-46 
& Ag 
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and WwW as 


Wy | 11-47 


The entire set of equations then takes the form 


wee 


A w = \p/ 11-48 


oe NUMERICAL SOLUTION OF THE MULTIGROUP EQUATIONS 
The basic equation II-48 can be solved by various direct or 


iterative methods. Direct methods generally are based on either 


eB) 
8, 
bh 


Gaussian elimination with pivoting, or on triangular decompositic 
the matrix of coefficients. Iterative methods include the Jacobi 
method, unextrapolated Liebmann or Gauss-Seidel method, various 
forms of the extrapolated Liebmann or successive orerrelaxation 
method and others. 

In this study, up to 25,000 algebraic equations were expected as 
a result of a sixteen group problem in two dimensions with 1600 mesh 
points. For this large number of equations, iterative methods are 
generally used rather than direct methods [8, 9]. Smith [9] states 
that one reason for this is that the iterative methods are more efficient 
in that they take advantage of the sparse matrix of coefficients with 
its numerous zero elements. 

In this study, the point successive overrelaxation method was 


used (See Appendix D). 


aS 





First the reactor materials and energy group structure are chosen. 
Next the cross section data is computed or as was the case in this 
study, taken from some reliable source {7 ] and used to compute the 
appropriate coefficients. Then initial estimates of the group fluxes 
designated Dy at all mesh points in the reactor except those lying 
on the extrapolated boundaries are made. These values are then used 
along with their respective coefficients to compute the initial source 
vectors Wg . the flux in group one at all mesh points in the 


reactor is then corrected by solving the set of equations 
A, go = Wi? 
Dy =A. 11-49 


Using the point successive overrelaxation method, this is done 
by solving the first of Equations II-49 for the central mesh point, 
mer, , C2,2) and then solving the remaining equations point-by- 
point, row-hy-row, at all other points no on the extrapolated boundary. 
The solution is denoted ZB . The remaining group equations are then 


similarly solved yielding the vector ZB where 
oD! 
l 
B=] 
? IT-50 
“oa t 
D x 


which represents the first iteration. 


This vector is then rescaled such that 


to, 2 0.8 )= CE", oD : I-51 
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which indicates that the flux vector on the zeroth iteration integrated 
over all spatial dimensions and all energy groups is equal to the flux 
vector after the first iteration integrated over all spatial dimensions 
and all energy groups multiplied by a scale factor [10]. When the 
process is carried out, it can be seen that after a few iterations the 
scale factor will approach a constant value. The reason for this 
as described by Lamarsh [11] is that each iteration is actually 
equivalent to one cycle of the chain reaction. The zeroth iteration 
of the fluxes bp ° then can be thought of as a set of initial conditions 
or the flux state at time t= 0. Eventually the flux will approach the 
funadmental eigenfunction as the higher harmonics die out. The 
fundamental will then vary with the number of iterations (time) 
increasing or’ decreasing depending on whether the system is 
elecritical or critical. The quantity *% can clearly be considered 
a measure of the criticality of the system. If % is less than unity 
then the group iluxes throughout the system are decrcasing with cach 
iteration and the system is subcritical andif 5 is greater than 
unity then the fluxes are increasing and the system is supercritical. 
5 is in fact the inverse of the effective multipication factor Boge. 
The iteration process described above is referred to as the inner 
iteration o flux iteration. If a critical assembly is desired, then 
an outer iteration must be performed wherein adjustments are made to 
the composition or size of the reactor in order to make the effective 
multiplication factor converge to unity. | 
Convergence of the flux vectors on the inner or flux iterations 
is determined by taking the absolute value of the difference between 
the flux at all mesh points and in all energy groups and comparing it 


with some desired value €, , such that 


N-l | 
| D Canige ~S(E,r,z) er Meise 
AS 


m,r,Z 


ZT 








Convergence of the outer or criticality iteration is determined by 


subtracting the scale factor % from unity so that 


1.0 - % is €, ; 11-53 


where €, is an arbitrary convergence criterion. 
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III. DESCRIPTION OF COMPUTER PROGRAM 


A. INPUT DATA 

In this section the input data is read in via punched cards as 
illustrated in Appendix E and printed in the output for verification. 
The input data along with their symbolic program designators include 
the following: the energy dependent source SS (I); the maximum values 
of the initial flux estimate S(I): the energy intervals DELE(I); the number 
of atoms per cubic centimeter H(I) of each element or isotope; the 
microscopic cross sections for each element or isotope for inelastic 
scattering SIGIN(I,J,K), transport SIGTR(I,J), fission SIGF(I,J), 
elastic removal SIGER(I,J), and capture SIGC(I,J); the normalized 
neutron fission spectrum CHI(I); and the average number of neutrons per 


fission NU(I,J). 


Eee QUTER ITERATION 

If a critical composition is the desired result, the program returns 
to this section after the iteration has converged (See SAMPLE COMPUTER 
PROGRAM pagell9). The iteration on composition then begins by 
adding to the volume of the reactor an incremental volume of vranium 
235. An equal volume of uranium 238 is subtracted and new values for 
N29 H(1) and n28 H(2) are computed. On the zeroth outer iteration 


this section of the program is by-passed. 


oe COMPUTATION OF MATRIX OF COEFFICIENTS 

Input data is used in this section to compute the elements of the 
coefficient matrix, the scattering matrix, and the fission source matrix. 
In the first part of this section, the various microscopic cross sections 
for each reactor component are multiplied by their respective atomic 
abundances N*, The values for each particular reaction (for example, 


capture or inelastic scatter from g to g + 1) are then summed. This 
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is done for each energy group. The summations are then combined 
with other input data and the coefficients are computed for each 
group equation. 

The equation denoted by statement number 17 in the sample 
program computes the diffusion coefficients for each group based on 
the simple diffusion theory approximation. The absorption coefficients 
for each group are computed in the equation denoted by statement 
number 175 by summing the contributions from all the absorption 
or loss reactions. The equation denoted by statement number 18 
computes the fission and down-scattering coefficients from all 
higher energy groups. This is accomplished by starting at the 
highest energy group in which the reactions originate and computing 
its contribution to all lower energy groups at the end of the slowing- 
down process. The fission source coefficients for a particular group 
and all lower energy groups are computed by the equation given in 
statement number 21. The above coefficients are completaly defined in 


Section II and Appendix B. 


D. INITIAL FLUX ESTIMATE AND BOUNDARY CONDITIONS 

In the cylindrical geometry the flux shape is initially estimated 
as the product of the relative group magnitude factor S{I), a zero 
order Bessel function, and a cosine function. The factors S(I) were 
determined for this study as suggested by Clark and Hansen [10]. 

The flux on the exterior boundaries is set equal to zero. At the 
interior boundaries the flux at the first mesh points on either side of 
the longitudinal and mid-plane axes are set equal, indicating that the 
flux has its maximum value along these axes. At the end of this section 
the initial flux estimate is printed out for verification and comparison 
with the final flux values. 

The program is designed to by-pass this section on all outer 


iterations. 
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BK INITIAL FISSION SOURCE INTEGRATION 
The initial flux estimate is used in this section to compute the 
fission source vector for the zeroth inner iteration. The program 


performs the following triple integration: 


ap (° 2 ce 
| | | i[ Bce,r,21 MCE) Dv: CE) VEE e) }dedr dz I-1 
o YO Vo bi 


by dividing the energy spectrum into discrete intervals DELE(I) 
corresponding to the energy groups (for this study the upper boundary 
of the highest energy level was taken as 10 Mev for use in this calcu- 
lation) and dividing the area into finite intervals corresponding to the 


mesh spacing GS. A summation process is then carried out. 


| 


COMPUTATION OF SOURCE VECTORS 

The components of the vector W for each inner iteration are computed 
by summing the products of the fission source or the combined fission 
source plus down-scattering source coefficients and the scaled flux 
values from the previous iteration. On the zeroth iteration the initial 

flux estimate is used with a scale factor of unity. The resulting numbers 
SUMR(I,R,Z) and SUML(I,R,Z) are the space dependent components of 


each group source vector Wag. 


.. POBUITION OF THE DiFPUSION EQUATION 

This section consists of one outer ener group loop with several 
space loops within it. Each group equation beginning with the highest 
group is solved at all spatial points starting at the longitudinal axis r= 0 
and proceeding column by column to within one mesh point of the extra- 


polated boundaries r = R-l and z = Z-1. 
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The first space loop will solve the diffusion equation at the 
geometric center of the reactor, i.e., point (2,2) (See Appendix C, 
Figure C-1.) if an extraneous point source is located in this position. 
This particular configuration was used in this study; however, extrane- 
ous sources with other spatial distributions could also be used with 
the appropriate changes to this section of the program. 

The second space loop involves the solution of the diffusion 
equation at all points on the longitudinal axis r= 0 but not occupied 
by - extraneous source. A special finite difference approximation 
for points on this axis is necessary due to a singularity caused by the 
radius being equal to zero (See Appendix A, page 66). 

The last spatial loop solves the diffusion equation at all remaining 
mesh points with the exception of those located on the exterior 
boundaries. 

Within each of these spatial loops, the flux value at each mesh 
point from the previous iteration is unscaled (TPHI1) and compared 
with the flux value at each mesh point in the present iteration to 
determine convergence as defined in Section II, page 27. The 


largest absolute value of these differences, 


TEST2 = DABS (TPHI2 - TPHI1) , IlI~2 


is stored as the quantity SAVE and may be printed with the output as 
desired. The scaled flux values of the previous iteration SPHI1, 
SPHI2, SPHI3 are used in the successive point overrelaxation scheme 


for solving the group equations at each mesh point. 


H. NTH FISSION SOURCE INTEGRATION (UNSCALED FLUX) 
In this section the triple integration process, Equation III-1, is 


performed with the new but as yet unscaled flux values. 
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I. COMPUTE SCALE FACTOR AND RENORMALIZE FLUX 

Once the fission source integration has been computed for the 
N-1 iteration FLUXI(NN-1) using the scaled or renormalized flux 
values and for the Nth iteration FLUXI(NN) using unscaled flux values, 


the scale may be determined from 
SCALF = DSQRT(FLUXI(NN-1)/FLUXI(NN)) III-2Z 


The flux values at all mesh points arid in all energy groups are then 
renormalized by multiplying by this scale factor. The scale factor 
SCALF is the principal eigenvalue of the system and its reciprocal is 


the effective multiplication constant SKEFF as explained in Section II, 


page 


J. NTH FISSION SOURCE INTEGRATION (SCALED FLUX) 
At this point the integration process of Section III(E) is performed 


Vide 


using the rencrmalized flux. On the subsequent iteration thi 


i) 


become the quantity FLUXI(NN-1) in the numerator of Equation III-2. 


K. END OF INNER ITERATION 

This statement (1801) is the controlling statement for the inner 
iteration. It can make the program exit to either the output section 
or the outer iteration loop depending upon whether or not a critical 
reactor is desired. The controlling factor may be either a specific 
number of inner iterations if experience has shown that the convergence 
criterion as defined in Section II, page 27 will be satisfied or another 
possibility would be to stop the inner iteration once the convergence 
criterion is met regardless of the number of inner iterations. The 
second method might be preferable once it had been determined that 
convergence to a reasonable criterion, i.e., €., would not take an 


excessive amount of computer time. 
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Pee OQULPUT SECTION 

The output is very flexible and may be easily changed to fit a 
particular problem or set of problems. Typical output data includes the 
following: scale factor, effective multiplication constant, flux 
values at various mesh points and various inner or outer iterations, 
convergence criterion SAVE (€, ) and overrelaxation factor w. In 
addition, graphical output in the form of graphs of flux shape compared 
with the cosine function and zero order Bessel function at various 
locations may be plotted and the convergence criterion SAVE versus 


iterations. Many other possibilities exist. 


mee oLART NTH OUTER ITERATION AND TEST FOR CRITICALITY 
In this section the scale factor is checked out for criticality 
by comparison with an arbitrarily selected value ( .0005 was chosen 


feretais study.) such that 


11.0-SCALF| << .0005 III-3 


If the above condition is not satisfied, then another outer iteration 
is made. The number of outer iterations may also be arbitrarily 
limited in this section. Once this limit, for example, ten outer 
iterations, is reached, the program will stop even though criticality 


has not been reached. 
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IV. RESULTS AND PROGRAM VERIFICATION 


If in the paragraphs that follow modifications to the basic core 
described below have been made, these modifications will be so stated. 
The results described and tabulated in this section were based ona 
typical homogeneous fast reactor core composed of the following 


materials with their respective volume fractions: 


U5 
U ~ 30% 
Na ~ 33% 
re ~ 27% 
pee. 8% 


This composition was chosen to coincide approximately with the 
RAPSODIE reactor core in reference [12]. The core dimensions were 
36 centimeters in radius and 72 centimeters in height. Mesh spacing 
was two centimeters in both the radial and axial directions. 
energy groups were used with the cross section data taken from 


ANL 5800 [7]. 


A. DETERMINATION OF OPTIMUM RELAXATION FACTOR 

The experimental method used for determining Eon consisted of 
two basic steps. First, an-analytical estimation was obtained based 
on formulation developed by Franke! [13] and Starling [14]. 

Frankel [13] derived the following equations based on the 
solution to Laplace's equation in a bounded region R in rectangular 


coordinates 


Wsysic (CL Wel 
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where @& is the smaller root of 
Oa = eee IV-2 


and tis given by 


where p and g are the number of mesh points in the x and y directions. 
Smith [9 ] used this same method to solve the similar Poisson's equa- 
tion in rectangular coordinates in a bounded region R. Starling [14] 
has developed an equation for t applicable to Laplace's equation in 


cylindrical coordinates and given by 


B.837h 
t = cos ce oe) 5 [V-4 


where Z and R are the height and radius of the cylinder respectivelv 
and his the mesh spacing. 

The table below gives the results found for these analytical 
solutions for Wags using a symmetric cylindrical section 36 centimeters 


in radius and 36 centimeters in halfi-height 


METHOD FOR FINDING "t" ene REMARKS 
Eq. IV-4 1.7888 Z=72cm. (full height) 
Eq. IV-4 1.7259 Z= 36 cm. (halt heighs) 
Eq. IV-3 1.7294 p = q= 20 (no. mesh pts. + 1) 
Eq. IV-3 1,704 p= q— 18 (io. meaner.) 


The second step in the experimental process was to try these values 
in the actual program and graph the convergence criterion versus the 
number of inner iterations for these values of eat along with other 
reasonable estimates to determine the best one. Figure IV-1 shows the 
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data which resulted from this study. It can be seen that for the values 
of wae equal — .72, 1.73 and 1.74, a convergence criterion value 
of less than 10 §is satisfied within fifty iterations. Based on this 
data, it was concluded that Equation IV-4 with Z equal to one-half 

the height of the cylinder, due to symmetry, gave the best value of 
Sot for this program. The program was modified to compute Mane 
automatically based on Equation IV-4. 

The above results were reasonable in view of the fact that 
Equation IV-4 was developed for use with Laplace's equation in 
cylindrical coordinates. The multigroup neutron diffusion equation 
used in this study in cylindrical geometry is of a similar form. One 
of the major differences is that the dependent variable @ isa 


function of both space and energy. 


feet LUX SHAPES FOR A TYPICAL NONCRITICAL CORE 

For the bare reactor described at the beginning of this section, the 
initial flux shape was estimated to be the product of the relative group 
weighting factor S(I), a cosine function and a zero order ordinary 
Bessel function of the first kind. This corresponds to the solution of 
Laplace's equation in a finite cylinder. Figure IV-2 shows the final 
flux shape in the axial direction at the center of the reactor after 
260 inner iterations with a convergence criterion of less than ie 
compared with the cosine Shee at the same location at zero iterations. 
The magnitudes were normalized to unity so that the shapes could be 
easily compared. It was found that the deviations were in the order of 
1 to 15 per cent with the per cent deviation increasing from the center 
outward. All three energy groups showed exactly the same normalized 
shapes after an equal number of iterations. The magnitudes of the 
PROUD fluxes are unknown since they depend on the power level at 


which the reactor is being operated. 
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The same comparison has been made in Figure IV-3 between the 
flux in the radial direction and the assumed initial Bessel function. 

A very slight deviation in shapes was noted here also indicating the 
close similarity between the diffusion equation and Laplace's 
equation. 

Examination of the data revealed that the final flux shapes were 
in effect established completely to seven decimal places after some 60 
inner iterations at which time the convergence criterion SAVE was less 
than 10” 

The program was run using a constant as the initial flux shape 
estimate rather than the cosine and Bessel function. The resulting 
final shapes for all three groups were exactly the same as those in 
Figures IV-2 and IV-3 after 90 inner iterations. Convergence was 
somewhat slower using this initial estimate. For example, after 50 
iterations, SAVE was equal to 3.2 X 10° versus 3.7 X io 2 for an 


initial estimate based on the cosine, Bessel function product. 


C. PROGRAM TEST WITH EXTRANEOUS POINT SOURCE 

The program was AG with a simulated extraneous Pu-Be point 
neutron source at the geometric center of the reactor. Two core 
compositions were used in these tests. 

The first core was made up of the typical composition described 
at the Neemnning Grins Sections ne tux shiapes ineine Tadiawana 
axial directions for the three groups were plotted on the same graph 
with the Bessel functions and cosine respectively. Figures IV-4 and 
IV-5 show the results of these comparisons after 90 inner iterations. 
Group one showed the largest deviation with the most pronounced 
effect nearest the source. This result was due to the normalized source 
neutron contribution to group one. Eighty-seven per cent of the neutrons 
from the Pu-Be source are born with energies lying in group one 


(1.35 Mev —* c& )as defined in ANL S5800L 7 J. 
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The second core was a solid iron cylinder with the same 
dimensions as the typical fast reactor core. The Pu-Be source 
was located at the center of the cylinder. The radial flux shapes for 
groups one and two were plotted in Figure IV-6 and the axial flux 
shapes for these two groups were plotted in Figure IV-7. It was noted 
that for these two groups the flux decayed very rapidly as the distance 
from the source increased; therefore, the ordinates of Figures IV-6 
and IV-7 are logarithmic. Figures IV-8 and IV-9 show the radial and 
aal flux for group three. In these graphs the flux shapes are com- 
pared with the cosine and Bessel functions after 90 inner iterations. 

The three group flux spectrum at the center of the two assemblies 
are shown in Figure IV-10. A relatively flat spectrum was observed for 
the core containing the fissile material while the solid iron core showed 


large relative differences amongthe three group fluxes. 


Dy Gm TiCALITY TEST 

Testof the program up to this point included only the inner iteration. 
Criticality was achieved by two methods of outer iteration. 

mine tirst metnou Comsisrved Of Perorming SUFIciICn! Inmer iteraieme 
for convergence. Seventy iterations were used resulting inan €, of 
less than 10°” for a mesh size of two centimeters. An outer iteration 
was then carried out wherein .5% by volume of wo was added before 
Seach outer iteration as described in Section III-B. This cycle was then 
repeated until the criticality criterion ( [1.0 - SKEFF\|<.0005 ) 
was satisfied. The optimum overrelaxation factor remains the same 
for each outer iteration. 

Four different mesh spacings were used in this section of the 
Piegram analysis to determine the effect, if any, mesh size had on 


criticality and flux shapes. Table IV-1 contains the data resulting 


from this study. The central flux shapes in the radial and axial 


directions at criticality for all mesh sizes tried were nearly coincident 
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with the Bessel function and cosine respectively. Table I1V-2 shows a 
maximum deviation of .71 per cent using a mesh size of two centimeters. 
Criticality was reached at the same core composition, i.e., 4.5% we 
for mesh spacings of .75 centimeters, 1.5 centimeters and 2 centimeters. 
With a mesh size of four centimeters, criticality was not achieved 

at this composition. 

The second method for arriving at a critical assembly consisted 
of fixing the core composition and varying the core dimensions in the 
outer iteration. A new optimum overrelaxation factor was computed 
for each outer iteration since the geometric parameters that were used 
to compute ae were changing. Three compositions corresponding 
to 18.9%, 20% and 27% enrichment were used. The results of this 
test are plotted in Figure IV-11 and compared with the bare core criti- 
cality curve of the RAPSODIE reactor in reference [11]. Differences 
noted are believed to be due to the following: a) compositions are not 
exactly the same: b) different cross section data was provabiy used 
though this was speculative since reference [11] did not state the 
source of their cross section data: and c) different programs were used 
to solve the multigroup diffusion equations. PRODII program was used 
for the RAPSODIE reactor. The factor above having the most influence 
was probably the large volume of iron. The absorption cross section 
of iron is fairly large and as seen in Table IV-3 is occupying a volume 
which in RAPSODIE contained 4% molybdenum, an unspecified quantity 
of iron in the form of stainless steel, and possibly some void space. 
The absorption cross section of these components is likely to be less 
than pure iron. 

The comparison of these curves shows that while a much larger 
critical reactor,assuming the same enrichment, would result from 
using the program developed in this study, the shapes of the curves 


are very similar. The fact that as enrichment was increased, the core 
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size and hence the critical mass was decreased indicates a correct 
imond., the ilux shapes found using this outer iteration scheme were 
also nearly coincident with the cosine and Bessel function curves. 

Table IV-4 shows an example of the per cent deviations experienced 

for a critical core assembly 32 centimeters in radius and 124 centimeters 


in height. 


CE. TESTS WITH SIXTEEN GROUP DATA 

Attempts to modify the program to solve a sixteen group problem 
were not successful. Both under- and over-relaxation factors were 
tried but satisfactory convergence of the inner iteration could not 
be achieved after over thirty minutes of computer time. An under- 
relaxation factor of . 1 was used and after 310 iterations and 34 
minutes, the convergence criterion was only of the order 10. 
Overrelaxing with a factor of 1.726 yielded a highly oscillatory flux 
behavior in all energy groups. All group fluxes were initially estimated 
to have maximum value of ten. The magnitude of oscillations recorded 
in the output varied between plus or minus eighty for the higher energy 
groups to plus or minus thirty for the lower energy groups with no 
noticeable damping after 310 iterations and some 34 minutes of 


computer time. No further tests were conducted. 
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MESH POINT INITIAL FLUX EVE ell UD.4 PER CENT 
LOCATION GROUPS. ls 27 Se GRO BS a 2 Se DEVAien 
Ca Ee eee ye 














(2,2) 1.000000 1.000000 

(42) .977535 .977479 

(6, 2) .911646 .911554 

(8,2) .806757 . 806654 

(10, 2) .669888 .669786 
(enn) 2) -510080 .510013 
(Ay 2) . 337808 .337782 
(16,2) mlo4t2d . 164154 
(18, 2) 0.0 0.0 

ie. 2) 000000 1.000000 --- 
(2,4) 994868 .994805 . 006 
(2,6) .979525 .979403 012 
ig? 8) .95414] .953952 .020 
(Ze, 1.0) .918959 .918717 .026 
2 12) 84a gals . 874059 mia2 
(2,14) .820767 . 820438 .040 
i, 16) .758757 .758407 .046 
(2,18) .688969 . 688603 053 
(2,20)  GaiZe. 01 So L1742 058 
(222) .528959 Poe os .065 
i), 24) .440393 .440076 a7 
(2, 26) .347303 .347033 .078 
(2,28) . 250656 .250440 .086 
(2, 30) JUS ag . 151289 .089 
(2,32) .050646 .050600 .091 
Table IV-4 
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V. CONCLUSIONS 


This study was undertaken to develop an easy to apply computer 
code for numerically solving the multigroup diffusion equations. The 
program was written in the FORTRAN IV language and all work was done 
on an IBM 360/67 computer. Applicability to other computer systems, 
while not a consideration of this study, is assumed to be practical. 

The unique feature of this program is that it will be readily 
available to and easily used by students at the Naval Postgraduate 
School. It is hoped that this study will complement future local 
experimental studies in fast neutron diffusion and also enable the 
student researchers to conserve valuable time which would be 
expended attempting to verify their experimental data with a myriad 
of complicated and unfamiliar multigroup diffusion programs. 

The point successive overrelaxation method of solution was 
successful with the three group data on the basis of the convergence 
criterion chosen, i.e., 


Oo C2ie 2) i Oo Cae a lo 
L 


Soe re 
Convergence of the three group criticality problem with five outer 


iterations took five minutes and forty seconds and 126K:words of 
core storage using a mesh size of two centimeters. Criticality was 
achieved with smaller mesh sizes, i.e., .75 centimeters and 1.5 
centimeters, but computation times were considered excessive especially 
for .75 centimeter mesh size which took about 33 minutes of computer 
time and over twice as much core storage space. 

Tests showed that the code was easily adaptable to three group 
solutions for typical fast reactor cores and reactor material assemblies 


with extraneous point sources. Flux shapes and energy spectrums from 
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these tests should provide a sound basis for future experimental flux 
measuring studies in fast reactommmac amar 

The close coincidence of the critical flux shapes with the cosine 
and Bessel functions indicates that for a three group problem if critical 
core size and composition are known, then the space-independent solu- 
tion using the buckling may be used. It was found, however, that 
this was not true for subcritical assemblies, i.e., the space and 
energy dependence of the flux are not separable. 

Further work with larger numbers of energy groups will have to be 
carried out to establish the program's ability to handle successfully many 
group calculations. The unsuccessful attempt at using sixteen groups 
was due to the failure to attain suitable convergence in a reasonable 
amount of time. A better technique than the point successive over- 
relaxation method is obviously needed when dealing with this number 


' of energy groups in order to improve the convergence time. 
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VI. RECOMMENDATIONS FOR FUTURE WORK 


In its present form the program has proven to be successful for 
the purpose for which it was intended. Future researchers might choose 


to expand or modify this program in one or more of the following areas. 


A, GEOMETRY AND SPATIAL DIMENSIONS 

Although the program is written for the cylindrical geometry in two 
Peatial dimensions (flux g is nota function of azimuth angle © ), 
it need not be so limited. The program could easily be adapted to the 
three dimensional cylindrical geometry as well as the two or three 
dimensional rectangular geometry. The addition of a third dimension in 
space would, of course, require considerably more computer core storage 
and hence the future user might have to limit the physical size and 
number of energy groups in a particular problem. 

Another possible additional feature might be to increase the 
number of regions by one or more to include blanket, reflector and 


shielding regions. 


B, CROSS SECTION DATA 

All the results for this program were obtained using three group 
ANL data. Conceivably future users would want to use the most 
current and widely accepted cross section available to them. For 
this reason the program could be modified to handle multigroup data of 
up to twenty-six groups or higher. Once again the user must bear in 
mind that as the number of energy groups is increased, core storage and 


computation time will also increase and probably not linearly. 


eo. INNER ITERATION METHODS 
The point successive overrelaxation iterative method gave 


satisfactory convergence rates for the three group problems used in 
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testing the program. However, more rapid convergence of the inner 
iteration is possible [15] using Hine ratteminan pompmenccessive 
overrelaxation methods and cyclic Chebyshev semi-iterative methods. 
If the future user desires to use more than three energy groups, he will 
almost certainly have to change the program to fit one of these methods 


or solution. 


1), TIME DEPENDENCE 

On a somewhat more far reaching level, this program might be 
expanded to include the time domain for the purpose of studying the 
Selamics Of a fast reactor. This could also include a study of space- 
dependent burnup and reactivity feedback effects. These topics are 
currently under investigation by those working on the development of 


the fast breeder reactor. 
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APPENDIX A 
DERIVATION OF LAPLACIAN OPERATOR IN CYLINDRICAL 
COORDINATES AND FINITE DIFFERENCE APPROXIMATIONS 


In order to numerically solve the set of partial differential 
equations II-31, a finite difference technique was chosen. In 


ry 
cylindrical geometry the Laplacian or WY operator is approximated 


by transforming the Cartesian coordinates to polar coordinates. 





PCx,y,Z2) =P(r,z,0) 


GEN ESE) Vo, © Singe 


r= (x?+y?) 9 @= tan! ¥ 
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Taking the partial derivations of rand © with respect to x and y in 


Equation A-1 results in 





- Te Zax ee a 
ax > 2 CXS y 2/2 —_> “Se cos © ) A-2 
a oY 4 ae ; 
8 _ Rea _ _ yy Aas = = peasy 
ae 2) Se ee r > a 
ee ee ee OES 

ae 7 Oe ee a DoS 

Now consider a function $¢+,2,@)  , whereinr and © are 


. functions of x and y through Equation A-1. Taking the first partial 
derivative of @ with respect to x, y, and z utilizing Equations A-2 
through A-5 gives 


bd, = Drvx + Bo Ox = Gr Coso-gdo SOO | A-6 


py = Dy as Deo Oy = Dr SiIn® + Be Sees re 


Dz = Dz Lo Bz : A-8 


63 





The second partial derivations are now obtained as follows through 


Equations A-6 through A-8: 
S Siné Sine 
ux = ( St COS © =~ 56 pe eee cin 


Sy. . 9 
Sir © Sint 6 
= Drv cos*6e +- Dr v =a Bee reo ya Dre eS 








Sine COsSe Sine cose 
r + 2 Po ae ) 


cela © GO= 2 
Pyy = Car Sino+ 2S cers) ( Dr. SinO +o ee) 


38 
g*6 Cos*e = 
= Derr Sint @ + Dr a ae Dee —) ae A-10 
COS6 SiN®e Cose Sin® 
+2prqo O88 SMe _ pg, L052 SS 
A-1li 


22 = a (Sz) = Dz2z 


Adding Equations A~9, A~-10 and A~-11, the Laplacian in cylindrical 


coordinates becomes 
V* B= Srr +t Oe +7 Boo + Pzz mee 


The partial derivatives are replaced by their respective central difference 
2 2 . | 
operators 16 with errors of order h where h is the square of the 
distance between adjacent mesh points (See Figure A-1.) resulting 
in 
7D = ee 3 B (r+, PA) 2G GAG) ) +Bcr-1,z,6) % 
\ | 
i ; ) Crtl,z,@) - G(r-l,z 0) 3 +o AEs 5 ® Cr, z,0-1) 
A-13 
| 
-2 Dor,z,e)+ Aly, z, O+1) 3 oe os ; D (r,Z+1,6) ~2B(r,Z,0) 


+ H(r,2-1, 8) 2 
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Pigure A=] 
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where the symbols used in Equation A-13 are shown in Figure A-l. 

Salvadori and Baron [16], as well as any other good text on 
numerical methods, derives finite difference operators with higher 
order errors, i.e€., more accurate approximations. 

It was found that this operator cannot be used fora mesh point 
lying on the longitudinal axis of the cylinder. At these points 
Equation A-12 has singularities in the second and third terms. In 
this study circular symmetry was assumed, thus Equation A-12 reduced 


to 
VoD = Dre + Br oD 27 A-14 


iinere still existed a singularity imthe second term; however, the 
problem under investigation also possessed symmetry with respect to 
the origin, i.e., @- =Oatr=0. Smith [9] shows that by 


Maclaurin's expansion 


Br Cr) = @rfo)+r Drr (O)+ar* Derr (O)+--- ) A-15 
however, D r(0) = 0, Therefore, as r approaches zero in the limit 
tne term “ Dr in Equation A-14 is the second partial derivative 
of @ atr=0 plus higher order terms. Equation A-14 can be written 


for mesh points on the longitudinal axis 


or in finite difference operator form 


wae D = =, $ Bor+i,z) ~ 26 Cr,z)+ Plr-1,z)3 
A-17 


+73 5 GUC 2+) ) eon, 2) CG ,z-1) 3 
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APPENDIX B 
DEFINITION OF COEFFMIG@ ic 


In order to simplify the writing of Equations II-35, they are first 
expanded in order to utilize the sixteen group data available in reference 
, | 


A typical equation, for example, group three, is written as follows: 


[ >. N° Conn j— 3 +X%3 00 ore,,) 4 D Caz) 
J={ 
ee a 





A3(\) 


e 
S 
ee ee ee 
Az(2) 


\ \ Cc , Pe ee eee ) 
S Ne) 
2 CYss 
es 
A3(3) 


\ 
+ Gen § 3 Craiay- dace, 293+-te Eogtrzen -2 456702) 


4 


G 

} : a 3 3 
Bs(r,2-1)3 | at N CTE, 3 a Jin, 32-3 + Oi1n > s—* 2 

>= 





J 3 3 3 
7 Tin, R—+» 5 sas in) ieee Se Sn alee a Toes ae 





A3(A) 


J 3 
i oahe a re B3Cv,Zz) 
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ues NCI Ds ae D3 (r,Z) se 56 
ee eee 


Cc A3 (5) 
Boa 
af j ya nN" (Ks Vie gine )| Die(r,Z)= Cn 


A3(18) 
This equation is then written in the form 


[ A3(4)+ A303) (42) da (rz) +[As@ (+ an) | 


B3(r+ijz> +[A33) (p32) ] B3Cz1) + e+ = A304, 042) 


B~-2 
x Bae) ype Ginza) - Aa (5) D3 Cpe) +- 5 2 
IAS CNS) Cee, 4) tr SAMs 2) 


The general form of this equation is 


[Ag (at) + Ag (g) (44)| Saenz) +[Ag(q) Rt an) 
A \ 

Go (r+1,Z) + L Aq Cg) C 2 soe | Pai, 2) s L Ag Cg ) 

(72) | fg Cr, z+.) + | Ag cg) (73) | Ba lr,z~i) B-3 

=AgM) DiCr,Zz)+ -+- + Ag Cg-1) Pag-i (+,Z) 

pera (972) Oo Cv ,2) 4-7 - + Ag Ce) Die Cr, 2) 


=p SCG 4 0 e 
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This general equation is condensed further to the form 


Al ggCr-,2)+A2 Gg lr,z) +A3 Gg(rsi,z) +A4 Balr, Z-1) 
ae AS Da Cr, Z+\) = AS, B,C, Z) mee +tASog- Do-t Cues) pan 


+ AFa Da OL eee Oc + AY ae Die gees, site DE Cs aD : 


In energy group one, Equation B-3 does not apply. The reason for this 
is that the terms on the right side of Equation B-3 are not fully applica- 
ble to group one. These terms have the following meaning: the sources 
of all neutrons that enter group g due to down scattering from all higher 
energy groups plus neutrons born into group g due to fissions in all 
higher energy groups are represented by the first set of terms on the 
right hand side of Equation B-3 [ Ag( | ) D | (r,z) + -- 
+Ag (9-1) Ba- (r,z) | ; 

In group one these terms must be zero since there are no higher energy 
groups. 

These apparent inconsistencies were accounted for in the computer 


program, i.e., these terms were set to zero in the first group equation. 
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APPENDIX C 
FINITE DIFFERENCE MESH POINT Wear) PaO hie aoe 


In order to take advantage of symmetry in the study and hence 
reduce the computer core storage space needed, the following mesh 
point scheme was adopted: a) all points GQ (1,x)= PD (3,x) and 

D foi) = D (x,3) due to symmetry, b) all extrapolated boundary 
points DB (L,x) = D (M,x) = 0 and c) no mesh points have zero 
index due to computer language incompatability. See Figure C-1 for 


further details. 
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SELES 
GAUSS-SEIDEL ITERATIVE METHOD 


In this technique, which is also called the unextrapolated 
Liebmann method, the point single step iterative method, or the method 
of successive displacements, the latest iterative values are used as 
soon as they become available. 


Consider the solution of Poisson's equation 


V°SZ =F xy) , : D-1 


in a rectangular region as in Figure D-1. Assume the solution is proceed- 


ing column by column from left to right. 


ia) 4 ‘o. jt2 Was e (irs. 


I+ 4 eae tee | 


- Py tt 
i+ | $ 
- PP | 
i ‘aly - an 


Be 


Figure D-1 
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If the point at which the equation is presently being solved is 
D (i,j), the iteration formula is 


) 


Br Ciyy = at OC, 3) + DW” Cit, §) 


+ ei ae) + DC, je eh? FCi,§) ; 


where the standard five-point central difference approximation is 


2. 
substituted for the V operator. 


POINT SUCCESSIVE OVERRELAXATION 


This method, also called the extrapolated Leibmann method, can 
be thought of as the sum of the Gauss-Seidel iterate Equation D-2 and 


the nth iterate. It is written as 
a Ci,i) = what Bo cinty) “+ b™ Cita, j) 
4 BCU yd OCH, et) + HFS) 2 | D-3 
+ (i-w) B” Ci,}) 


where w is the overrelaxation or acceleration factor. 
The optimum value of this factor in generalis \<Wopt<2 and 


can be computed from 


a ee D-4 


opt i A l~ 2 
where A, is the largest eigenvalue of the problem resulting from a 
solution using the method of successive displacements, i.e., the 
Jacobian method. In practice, however, the eigenvalues are seldom 
known and finding them by solving the problem by successive displace- 
ments can be very difficult due to slow convergence. References 8, 9, 
Oy 17 “eooniene an experimental process for determining Wert for a 
particular problem. In this process various values of w are tried until 


one is found that gives the most rapid convergence. This is then ceic:mined 


ee: 





to be eee Reference 16 suggests two other methods for determining 
a 

The experimental approach described above was used in this 
study. | 

The rapid convergence experienced with this method on many 
different problems has proved far superior to the ordinary Gauss-Seidel 
or successive displacements method without overrelaxation. References 
8, 17 contain some examples. 


References 8, 9, 10, 17 relate more explicit details on the 


method of successive overrelaxation in varying degrees. 
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APPENDIX _E 
INPUT DATA SYMBOLS_ 
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